In [1]:
%pylab inline
import matplotlib.pyplot as plt
#from IPython.core.display import HTML
Tight bindig model of diamond and zincblende type semiconductors¶
In [2]:
def matafn(h):
'Array flatten a matrix list of appropriate dimensions'
H=hstack(h[0])
for sor in range(1,len(h)):
H=vstack([H,hstack(h[sor])])
return H
In [3]:
# Az abra kimentesehez az alabbiakat a plt.show() ele kell tenni!!!
#savefig('fig_rainbow_p2_1ray.pdf'); # Abra kimentese
#savefig('fig_rainbow_p2_1ray.eps'); # Abra kimentese
# Abra es fontmeretek
xfig_meret= 9 # 12 nagy abrahoz
yfig_meret= 6 # 12 nagy abrahoz
xyticks_meret= 15 # 20 nagy abrahoz
xylabel_meret= 20 # 30 nagy abrahoz
legend_meret= 21 # 30 nagy abrahoz
In [4]:
One=matrix([[ 1, 0],[0, 1]])
Zero = matrix([[ 0, 0],[0, 0]])
Ham1 = zeros(shape=(8,8), dtype=complex)
Parameters (from Yu & Cardona):¶
\begin{eqnarray} V_{ss} &=& 4 V_{ss\sigma}, \hspace{5mm} V_{sp} = \frac{4}{\sqrt{3}}\, V_{sp\sigma}, \nonumber \\ V_{xx} &=& \frac{4}{3}\, ( V_{pp\sigma}+ 2 V_{pp\pi}), \hspace{5mm} V_{xy} &=& \frac{4}{3}\, ( V_{pp\sigma} - V_{pp\pi}). \end{eqnarray}In [5]:
def Ham_op(k,params):
Ham = zeros(shape=(8,8), dtype=complex)
# notation from Yu& Cardona:
EsA, EsB, EpA, EpB, Vss, Vsp, Vxx, Vxy = params
d1 = 1/4*array([1,1,1])
d2 = 1/4*array([1,-1,-1])
d3 = 1/4*array([-1,1,-1])
d4 = 1/4*array([-1,-1,1])
expk1=exp(1j*dot(k,d1))
expk2=exp(1j*dot(k,d2))
expk3=exp(1j*dot(k,d3))
expk4=exp(1j*dot(k,d4))
# notation from Yu& Cardona:
g1 = (expk1 + expk2 + expk3 + expk4)/4
g2 = (expk1 + expk2 - expk3 - expk4)/4
g3 = (expk1 - expk2 + expk3 - expk4)/4
g4 = (expk1 - expk2 - expk3 + expk4)/4
# Basis of Ham is based on
# Rolf Enderlein, Norman J. Horing: Fundamental of Semiconductor Physics and Devices
# (s1, x1, y1, z1, s2, x2, y2, z2)
### 2nd nearest hopping can be included:
#see /home/cserti/okt/Szilfiz/Papers_books/Tight_binding_model/TB_Chadi_Cohen_Phys_Satus_Solidi_1975.pdf
Uxx=0 #Uxx=-1.46 now it is swithed off
gg= cos(k[1])*cos(k[2])
h11 = matrix([[EsA+Uxx*gg,0,0,0],
[0,EpA,0,0],
[0,0,EpA,0],
[0,0,0,EpA]])
h22 = matrix([[EsB,0,0,0],
[0,EpB+Uxx*gg,0,0],
[0,0,EpB,0],
[0,0,0,EpB]])
h12 = matrix([[Vss*g1, Vsp*g2, Vsp*g3, Vsp*g4],
[-Vsp*g2, Vxx*g1, Vxy*g4, Vxy*g3],
[-Vsp*g3, Vxy*g4, Vxx*g1, Vxy*g2],
[-Vsp*g4, Vxy*g3, Vxy*g2, Vxx*g1]])
h21 = h12.H
Ham = matafn([[h11,h12],[h21,h22]])
return Ham
In [6]:
# parameters for Silicon
# energies in units of eV
# For parameters see Table 2.26 in Yu & Cardona book
EsA= -13.55;
EsB = EsA;
EpA = -6.52+0.17; # -6.52
EpB = EpA
Vss_YC = -8.13;
Vsp_YC = 5.88;
Vxx_YC = 3.17; ## itt lehet, hogy eliras van,
## mert a Chadi and M. L. Cohen cikkre hivatkoznak
Vxy_YC = 7.51; ### itt mas cikkben van egy -1 szorzo
params_Si_YC = [EsA, EsB, EpA, EpB, Vss_YC,Vsp_YC,Vxx_YC,Vxy_YC]
# For parameters see Table 1. in file 'TB_Chadi_Cohen_Phys_Satus_Solidi_1975.pdf'
# D. J. Chadi and M. L. Cohen, "Tight-binding calculations of the
# valence bands of diamond and zincblende crystals,"
# Phys. Status Solidi (b) 68, 405 (1975).
EsA= -13.55;
EsB = EsA;
EpA = -6.52+0.17; # -6.52
EpB = EpA
Vss_CC = -8.13;
Vsp_CC = 5.88;
Vxx_CC = 1.71; ## itt lehet, hogy eliras van a Chadi and M. L. Cohen cikkben,
## mert nem jo eredmenyt ad a program.
Vxy_CC = 7.51; ### itt mas cikkben van egy -1 szorzo
params_Si_CC = [EsA, EsB, EpA, EpB, Vss_CC,Vsp_CC,Vxx_CC,Vxy_CC]
## For parameters see Table 2.7 in Enderlein & Horing book:
#Vss = 4*Vsssig
#Vsp = 4*Vspsig/sqrt(3)
#Vxx = 4/3*(Vppsig+2*Vpppi)
#Vxy = 4/3*(Vppsig-Vpppi)
Vsssig_EH = -1.93
Vspsig_EH = 2.54
Vppsig_EH = 4.47
Vpppi_EH = -1.12
params_Si_EH = [EsA, EsB, EpA, EpB, 4*Vsssig_EH,4*Vspsig_EH/sqrt(3),
4/3*(Vppsig_EH + 2*Vpppi_EH),4/3*(Vppsig_EH - Vpppi_EH )]
## For parameters see file 'TB_Si_cinquanta.pdf'
Vsssig_cin = -1.774
Vspsig_cin = -0.8473
Vppsig_cin = 1.531
Vpppi_cin = -0.7165
params_Si_cin = [EsA, EsB, EpA, EpB, 4*Vsssig_cin,4*Vspsig_cin/sqrt(3),
4/3*(Vppsig_cin+2*Vpppi_cin),4/3*(Vppsig_cin-Vpppi_cin)]
print("Yu & Cardona: \n")
print("Ep-Es = ",EpA-EsA)
print("EsA, EsB, EpA, EpB, Vss,Vsp,Vxx,Vxy =",params_Si_YC)
print("\nChadi & Cohen: \n")
print("EsA, EsB, EpA, EpB, Vss,Vsp,Vxx,Vxy =",params_Si_CC)
print("\nEnderlein & Horing: \n")
print("EsA, EsB, EpA, EpB, Vss,Vsp,Vxx,Vxy =",params_Si_EH)
print("\nfrom TB_Si_cinquanta.pdf: \n")
print("EsA, EsB, EpA, EpB, Vss,Vsp,Vxx,Vxy =",params_Si_cin)
In [7]:
## fcc lattice:
# Silicon, params_Si
params = params_Si_YC # from Yu & Cardona
#params = params_Si_CC # from Chadi & Cohen
#params = params_Si_Enderlin # from Enderlein & Horing
#params = params_Si_cinquanta # from TB_Si_cinquanta.pdf
Npoints = 20; ## hany pontbol keszul a gorbe
ymax = 24.5;
## unit cell vectors of the reciprocal vectors:
## in units of 2 pi/a
b1 = 2*pi*array([-1,1,1])
b2 = 2*pi*array([1,-1,1])
b3 = 2*pi*array([1,1,-1])
## high symmetry points in the Brillouin zone
## FONTOS: a szomszedos spec. k pontokat kell venni
## (ezek most nem a fenti abranak megfelelo pontok, de szomszedos pontok,
## es igy jo az abra)
Gp = 2*pi*array([0,0,0])
Xp = 2*pi*array([0,1,0])
Wp = 2*pi*array([1/2,1,0])
Kp = 2*pi*array([3/4,3/4,0])
Lp = 2*pi*array([1/2,1/2,1/2])
Up = 2*pi*array([1/4,1,1/4])
## G-X-W-L-G-K-X line,,,
#specpoints = array([Gp,Xp,Wp,Lp,Gp,Kp,Xp])
# L-G-X-U-G line,,,
#specpoints = array([Lp,Gp,Xp,Up,Kp,Gp]) # U,K -S - X
#specpoints = array([Lp,Gp,Xp,Up,Gp]) # U,K -S - X
specpoints = array([Lp,Gp,Xp,Up,Kp,Gp])
## FIGYELEM: fcc racs eseten az irodalomban az U,K jeloles azt jelenti, hogy
## U-K vonal menten NEM rajzoljuk ki a savot!!!!!
### kiiratashoz, az x-tengelyen:
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Wp):"$W$",
tuple(Lp):"$L$",tuple(Up):"$U,K$",tuple(Kp):"$\Gamma$"}
klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
if not i==3: ########## Az U,K pontok egybe!!!!!
sl = specpoints[i+1]-specpoints[i]
sl_hossz = sqrt(dot(sl,sl))
dk = sl/Npoints
for num in arange(0,Npoints):
# k vector along the line given by speclines:
klist.append(specpoints[i] + num*dk)
# length of the k vector and shifted:
tmp += sqrt(dot(dk,dk))
kk.append(tmp)
specpos.append(tmp)
klist.append(klist[-1]+dk)
E_shift = eigvalsh(Ham_op(Gp, params))[1]
enlist = []
for kv in klist:
enlist.append(eigvalsh(Ham_op(kv, params))-E_shift)
#enlist.append(eigvalsh(Ham_op(kv, params_Si_Enderlin))-E_shift)
enlist= array(enlist)
# drawing the figure
figsize(8,8)
subplot(111)
# len(enlist[0,:])
for ii in range(0,8):
plot(kk,enlist[:,ii],'r')
E_x_axis = (eigvalsh(Ham_op(Gp, params))[0]-E_shift)*1.2
i = 0
for xc in specpos:
axvline(x=xc,color='b')
text(xc,E_x_axis,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
i=i+1
ylabel(r'$E(\mathbf{k}) \,\, (eV)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
yticks(arange(-12,12,2))
xlim(0,specpos[-1])
ylim(-13,12)
cim = "Si, diamond, TB approx."
title(cim,fontsize=xylabel_meret)
grid();
#savefig('Fig_Si_TB_LGXU-KG_from_Yu_Cardona.png'); # Abra kimentese
Check, at the $\Gamma$ point the eigenenergies $E_i(\mathbf{0}), i=1,2, \dots 8$ can be obtained analytically (Enderlein & Horing):¶
$E_1(\mathbf{0}) = E_s+V_{ss}$
$E_{2,3,4}(\mathbf{0}) = E_p-V_{xx}$
$E_{5,6,7}(\mathbf{0}) = E_p+V_{xx}$
$E_8(\mathbf{0}) = E_s-V_{ss} $
In [8]:
## a sorrend nem ok, de a szamok igen.
E_at_G_YC = [EsA+Vss_YC,EpA-Vxx_YC,EpA-Vxx_YC,EpA-Vxx_YC,
EpA+Vxx_YC,EpA+Vxx_YC,EpA+Vxx_YC,EsA-Vss_YC]
for i in range(8):
print("i = %d, E_i = %8.3f, from code = %8.3f " % (i+1, E_at_G_YC[i] - E_shift,
eigvalsh(Ham_op(Gp, params))[i]- E_shift))
In [9]:
eigvalsh(Ham_op(Gp, params)), eigvalsh(Ham_op(Kp, params))
Out[9]:
In [10]:
eigvalsh(Ham_op(Gp, params))[1]-eigvalsh(Ham_op(Gp, params))[0]
Out[10]:
In [11]:
eigvalsh(Ham_op(Gp, params))[4]-eigvalsh(Ham_op(Gp, params))[0]
Out[11]:
In [12]:
eigvalsh(Ham_op(Gp, params))[5]-eigvalsh(Ham_op(Gp, params))[1]
Out[12]:
In [13]:
params
Out[13]:
In [14]:
eigvalsh(Ham_op(Gp, params))
Out[14]:
In [ ]: